P03 - Study of the energy resolution and uncertainties of photon detectors using Bayesian methods.

P03 - Study of the energy resolution and uncertainties of photons in Ge detector using Bayesian methods.

Germanium detectors have wide fields of application for γ- and X-ray spectrometry thanks to their excellent energy resolution. The energy resolution of these detectors is defined as the width of the detected energy spectra peaks (FWHM); it depends on

• the statistics of the charge creation process

• the properties of the detector, and primarily its charge collection efficiency

• the electronics noise

The resolution can be expressed as the squared sum of two terms

\[FWHM = \sqrt{w_d^2 + w_e^2}\]

where the first term depends on the detector properties as

\[w_d = 2 \sqrt{2 \ln{2} · F · Eγ · w }\]

with F the Fano factor1 , Eγ the energy of the photon deposited energy and w is the electron-hole production energy threshold in germanium (w ∼ 3 eV)[1] The other term in eq. 1, we is connected with the readout electronics and depends on the detector capacitance, the size of the detector and the bias voltage. The following plot shows an uncalibrated energy spectrum collected with a Germanium detector irradiated by a combination of three sources: 241 Am, 60 Co and 137 Cs[2].

According to [2], the source nuclides emit the following photons:

Nuclide	Photon energy (keV)
241 Am	59.54
137 Cs	661.66
60 Co	1173.24, 1332.51

and these are the first four peaks (starting from the left side) visible in the figure. Similar spectra have been collected with other gamma sources (i.e. Th-228).

My Image

1. Infere the FWHM and the Centroids of each peak

Using statistical methods similar to that presented during the course, infere the FWHM of each γ peak for all available γ sources

First of all, we import the data and plot the spectra of the three sources.

library(rjags)

#extract dataset from file and create a dataframe with them
counts<-read.table("IC48A_AmCsCo.txt", skip =1,header = TRUE)
counts<-na.omit(counts)
dataset <- data.frame(counts)
colnames(dataset) <- c("ADC_channel", "Counts")
dataset<-dataset[(dataset$Counts>=0),]

# Set the file path and name for the output PNG file
output_file <- "plot_Ge.png"

# Open the PNG device
png(file = output_file, width = 800, height = 600)  # Adjust width and height as needed

#plot the dataset
plot(dataset$ADC_channel,log(dataset$Counts), type = 's', col='red', main = 'Am−Cs−Co spectra', xlab = 'ADC channel [#]', ylab = "Log Counts [#]",cex.main = 1.5)

# Close the PNG device and save the plot as a PNG file
dev.off()

null device 
          1 

The parameters of the peaks are inferred with a JAGS model.

Generative model:

Signal: Gaussian fit

\[\text{Signal}(x |\text{N},\mu,\sigma) = \frac{N}{\sqrt{2 \pi}\sigma} e^{[ -\frac{(x-\mu)^2}{2\sigma^2}]} \] - \(x\): ADC Channel

• \(N\): amplitude of the signal

• \(\mu\): centroid

• \(\sigma\): standard deviation

Background: Sigmoid fit

\[\text{Bkg_sigmoid}(x |c_{1},c_{2},\mu_{\sigma},\text{offset}) = \frac{c_{1}}{1+ e^{c_{2}(x-\mu_{\sigma})} } + \text{offset}\]

• \(c_{1}\): amplitude coefficient

• \(c_{2}\): smoothness coefficient

• \(\mu_{\sigma}\): centre

Background: Linear Exponential

\[\text{Bkg_exp}(x |\alpha,\text{offset}) = e^{\tan[\alpha]x + \text{offset}} \]

• \(\alpha\): slope [rad.]

• \(offset\): offset

Measurement model: Poissonian distribution

\[\text{Poiss}(n) = \frac{\lambda^{n}e^{-\lambda}}{n!} \]

• \(\lambda\) = value of the generative model (signal + background)

• n = measured value (counts on the ADC Channel)

As Priors distribution for the Bayesian Inferesence we use exclusively Uniform distribution. We change the domain range for every peaks.

In the first peak we use as background a sigmoid function:

\[\text{Bkg_sigmoid}(x |c_{1},c_{2},\mu_{\sigma},\text{offset}) = \frac{c_{1}}{1+ e^{c_{2}(\mu_{\sigma} - x)} } + \text{offset}\]

241-Am En = 59.54 keV

# FIRST PEAK

#select the first peak
peak1 = dataset[(dataset$ADC_channel>=195) & (dataset$ADC_channel<=250),]
#max1=peak1[which.max(range1$Counts),"ADC_channel"]

#plot it
plot(peak1$ADC_channel,log(peak1$Counts),  col='red', main = 'Am−Cs−Co spectra', xlab = 'ADC channel [#]', ylab = " Log Counts [#]",cex.main = 1.5,)

x <- peak1$ADC_channel

#define the first model with a sigmoid background
exp.background.model.1 <- "model{
    #data likelihood
    for (i in 1:length(x)){

         y[i] ~ dpois(points[i])
         
         points[i] <- N/sqrt(2*3.14*sigma*sigma)*exp((-(x[i]-mu)*(x[i]-mu))/(2*sigma*sigma)) - (c.1/(1+exp(c.2*(x[i]-c.mu)))) + offset
                                
    }

    #priors for the generative model

    N ~ dunif(10000,80000)           
    mu ~ dunif(200,220)        
    sigma ~ dunif(0.5, 3)            
    c.1 ~ dunif(100, 300)          
    c.2 ~ dunif(0,1)      
    c.mu ~ dunif(180,240)            
    offset ~ dunif(100,300)         

}"

peak1.jm <- jags.model(textConnection(exp.background.model.1), 
    data = list(x = peak1$ADC_channel, y = peak1$Counts), n.chain = 3, quiet = TRUE)

update(peak1.jm , 1000, progress.bar = "none")

peak1.cs <- coda.samples(model = peak1.jm, variable.names = c("N", "mu", "sigma", "c.1", "c.2", "c.mu", "offset"), n.iter = 10000, thin = 15, progress.bar = "none") 

summary(peak1.cs)

Iterations = 2015:11990
Thinning interval = 15 
Number of chains = 3 
Sample size per chain = 666 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

            Mean        SD  Naive SE Time-series SE
N      1.854e+04 1.409e+02 3.1531487      3.3901930
c.1    2.410e+02 3.716e+00 0.0831370      1.1251305
c.2    2.963e-01 2.410e-02 0.0005391      0.0021578
c.mu   2.125e+02 2.981e-01 0.0066683      0.0334945
mu     2.062e+02 9.733e-03 0.0002178      0.0002179
offset 2.398e+02 3.435e+00 0.0768487      1.0579940
sigma  1.306e+00 7.815e-03 0.0001748      0.0001948

2. Quantiles for each variable:

            2.5%       25%       50%       75%     97.5%
N      1.827e+04 1.845e+04 1.854e+04 1.863e+04 1.882e+04
c.1    2.338e+02 2.380e+02 2.415e+02 2.442e+02 2.468e+02
c.2    2.523e-01 2.797e-01 2.957e-01 3.106e-01 3.473e-01
c.mu   2.120e+02 2.123e+02 2.125e+02 2.128e+02 2.131e+02
mu     2.061e+02 2.062e+02 2.062e+02 2.062e+02 2.062e+02
offset 2.332e+02 2.371e+02 2.403e+02 2.425e+02 2.452e+02
sigma  1.291e+00 1.301e+00 1.306e+00 1.311e+00 1.321e+00

# plot of the coda samples
plot(peak1.cs)

peak1_mcmc <- as.mcmc(do.call(rbind, peak1.cs))
#head(peak1_mcmc)

N <- mean(peak1_mcmc[,1])
c.1 <- mean(peak1_mcmc[,2])
c.2 <- mean(peak1_mcmc[,3])
c.mu <- mean(peak1_mcmc[,4])
mu <- mean(peak1_mcmc[,5])
offset <- mean(peak1_mcmc[,6])
sigma <- mean(peak1_mcmc[,7])


# fit
pred_peak1 <- N/sqrt(2*3.14*sigma*sigma)*exp((-(x-mu)*(x-mu))/(2*sigma*sigma)) - (c.1/(1+exp(c.2*(x-c.mu)))) + offset


plot(peak1$ADC_channel,peak1$Counts,  col='red', main = 'Am−Cs−Co spectra, Peak 1', xlab = 'ADC channel [#]', ylab = "Counts [#]",cex.main = 1.5)
lines(peak1$ADC_channel, pred_peak1)

For the rest of the peaks we use a Linear exponential background in the generative model:

\[\text{Bkg_exp}(x |\alpha,\text{offset}, \mu) = e^{\tan(\alpha)(x - \mu) + \text{offset}} \]

137-Cs En = 661.66 keV

#SECOND PEAK

# range of the second peak
peak2 = dataset[!(dataset$ADC_channel<=2200) & (dataset$ADC_channel<=2350),]
#max2=peak2[which.max(range2$Counts),"ADC_channel"]
#plot the second peak
plot(peak2$ADC_channel,log(peak2$Counts),  col='red', main = 'Am−Cs−Co spectra', xlab = 'ADC channel [#]', ylab = " Log Counts [#]",cex.main = 1.5)

# rjags model as .bugs file
# peak describe a gaussian
# with an exponential background exp(m(x - mu) + offset)
exp.background.model.2 <- "model{
    #data likelihood
    for (i in 1:length(x)){

         y[i] ~ dpois(points[i])
         
        points[i] <- N/sqrt(2*3.14*sigma*sigma)*exp((-(x[i]-mu)*(x[i]-mu))/(2*sigma*sigma)) +exp(tan(slope.alpha)*(x[i]-mu) + offset)  # line generative model
        
    }

    #priors for the generative model

    N ~ dunif(100000,200000)            
    mu ~ dunif(2270,2295)        
    sigma ~ dunif(0.01,3)             
    offset ~ dunif(0,100)        
    slope.alpha ~ dunif(-3.14/2, 0)  


}"

# initialization of the model
peak2.jm <- jags.model(textConnection(exp.background.model.2), 
    data = list(x = peak2$ADC_channel, y = peak2$Counts), n.chain = 3, quiet = TRUE)

# burn-in
update(peak2.jm , 1000, progress.bar = "none")

# creates the coda sample
peak2.cs <- coda.samples(model = peak2.jm, variable.names = c("N", "mu", "sigma", "offset", "slope.alpha"), n.iter = 10000, thin = 15, progress.bar="none") 

summary(peak2.cs)

Iterations = 2015:11990
Thinning interval = 15 
Number of chains = 3 
Sample size per chain = 666 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

                  Mean        SD  Naive SE Time-series SE
N            1.610e+05 4.013e+02 8.978e+00      8.977e+00
mu           2.279e+03 6.121e-03 1.369e-04      1.370e-04
offset       4.637e+00 8.931e-03 1.998e-04      1.997e-04
sigma        2.379e+00 4.624e-03 1.035e-04      1.079e-04
slope.alpha -4.165e-03 1.829e-04 4.091e-06      4.093e-06

2. Quantiles for each variable:

                  2.5%        25%        50%        75%      97.5%
N            1.602e+05  1.607e+05  1.610e+05  1.613e+05  1.618e+05
mu           2.279e+03  2.279e+03  2.279e+03  2.279e+03  2.279e+03
offset       4.619e+00  4.631e+00  4.637e+00  4.643e+00  4.654e+00
sigma        2.370e+00  2.376e+00  2.379e+00  2.383e+00  2.388e+00
slope.alpha -4.522e-03 -4.289e-03 -4.165e-03 -4.040e-03 -3.815e-03

# plot of the coda samples
plot(peak2.cs)

# merge the chains on a mcmc objects
peak2_mcmc <- as.mcmc(do.call(rbind, peak2.cs))
#head(peak2_mcmc)

# mean values of the parameters form the mcmc object
peak2_params <- c("N" = mean(peak2_mcmc[,1]),"mu" = mean(peak2_mcmc[,2]), "offset" = mean(peak2_mcmc[,3]), "sigma" = mean(peak2_mcmc[,4]), "alpha" = mean(peak2_mcmc[,5]))
#peak2_params[4] *  2.3

# fit of the gaussian peaks + the exponential background
pred_peak2 <-  peak2_params[1] /sqrt(2*3.14*peak2_params[4]*peak2_params[4])*exp((-(peak2$ADC_channel-peak2_params[2])*(peak2$ADC_channel-peak2_params[2]))/(2*peak2_params[4]*peak2_params[4])) +exp(tan(peak2_params[5])*(peak2$ADC_channel - peak2_params[2]) + peak2_params[3]) 

#plot of the data and the fit
plot(peak2$ADC_channel,peak2$Counts,  col='red', main = 'Am−Cs−Co spectra, Peak 2', xlab = 'ADC channel [#]', ylab = "Counts [#]",cex.main = 1.5)
lines(peak2$ADC_channel, pred_peak2)

NA
NA

60-Co En = 1173.24 keV

# THIRD PEAK

# range of the third peak
peak3 = dataset[!(dataset$ADC_channel<=3950) & (dataset$ADC_channel<=4100),]
#max3=peak3[which.max(range3$Counts),"ADC_channel"]

#plot the third peak
plot(peak3$ADC_channel,log(peak3$Counts), col='red', main = 'Am−Cs−Co spectra', xlab = 'ADC channel [#]', ylab = "Log Counts [#]",cex.main = 1.5)

# rjags model as .bugs file
# peak describe a gaussian
# with an exponential background exp(m(x - mu) + offset)
exp.background.model.3 <- "model{
    #data likelihood
    for (i in 1:length(x)){

         y[i] ~ dpois(points[i])
         
        points[i] <- N/sqrt(2*3.14*sigma*sigma)*exp((-(x[i]-mu)*(x[i]-mu))/(2*sigma*sigma)) +exp(tan(slope.alpha)*(x[i]-mu) + offset)  # line generative model
        
    }

    #priors for the generative model

    N ~ dunif(10000,90000)            
    mu ~ dunif(4025,4055)        
    sigma ~ dunif(0.1,5)             
    offset ~ dunif(0,400)        
    slope.alpha ~ dunif(-3.14/2, 0)  

}"

# initialization of the model
peak3.jm <- jags.model(textConnection(exp.background.model.3), 
    data = list(x = peak3$ADC_channel, y = peak3$Counts), n.chain = 3, quiet = TRUE)

# burn-in
update(peak3.jm , 1000, progress.bar="none")

# creates the coda sample
peak3.cs <- coda.samples(model = peak3.jm, variable.names = c("N", "mu", "sigma", "offset", "slope.alpha"), n.iter = 10000, thin = 15, progress.bar="none") 

summary(peak3.cs)

Iterations = 2015:11990
Thinning interval = 15 
Number of chains = 3 
Sample size per chain = 666 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

                  Mean        SD  Naive SE Time-series SE
N            6.136e+04 2.492e+02 5.574e+00      5.573e+00
mu           4.040e+03 1.388e-02 3.105e-04      3.104e-04
offset       3.957e+00 1.340e-02 2.998e-04      2.991e-04
sigma        3.323e+00 1.056e-02 2.363e-04      2.482e-04
slope.alpha -3.665e-03 2.594e-04 5.803e-06      5.806e-06

2. Quantiles for each variable:

                  2.5%        25%        50%        75%      97.5%
N            6.087e+04  6.120e+04  6.135e+04  6.152e+04  6.186e+04
mu           4.040e+03  4.040e+03  4.040e+03  4.040e+03  4.040e+03
offset       3.931e+00  3.948e+00  3.957e+00  3.966e+00  3.982e+00
sigma        3.301e+00  3.315e+00  3.323e+00  3.330e+00  3.343e+00
slope.alpha -4.169e-03 -3.840e-03 -3.660e-03 -3.491e-03 -3.163e-03

# plot of the coda samples
plot(peak3.cs)

# merge the chains on a mcmc objects
peak3_mcmc <- as.mcmc(do.call(rbind, peak3.cs))
#head(peak3_mcmc)

# mean values of the parameters form the mcmc object
peak3_params <- c("N" = mean(peak3_mcmc[,1]),"mu" = mean(peak3_mcmc[,2]), "offset" = mean(peak3_mcmc[,3]), "sigma" = mean(peak3_mcmc[,4]), "alpha" = mean(peak3_mcmc[,5]))

# fit of the gaussian peaks + the exponential background
pred_peak3 <-  peak3_params[1] /sqrt(2*3.14*peak3_params[4]*peak3_params[4])*exp((-(peak3$ADC_channel-peak3_params[2])*(peak3$ADC_channel-peak3_params[2]))/(2*peak3_params[4]*peak3_params[4])) +exp(tan(peak3_params[5])*(peak3$ADC_channel - peak3_params[2]) + peak3_params[3]) 

#plot the fit and the measurements
plot(peak3$ADC_channel, peak3$Counts,  col='red', main = 'Am−Cs−Co spectra, Peak 3', xlab = 'ADC channel [#]', ylab = "Counts [#]",cex.main = 1.5)
lines(peak3$ADC_channel, pred_peak3)

#peak3.acf <- acf(peak3_mcmc, plot = F) #plot = F will not show the plot of the ACF.

60-Co En = 1332.51 keV

# FOURTH PEAK 

# range of the fourth peak
peak4 = dataset[!(dataset$ADC_channel<=4500) & (dataset$ADC_channel<=4700),]
#max4=peak4[which.max(range4$Counts),"ADC_channel"]

#plot the third peak
plot(peak4$ADC_channel, log(peak4$Counts), col='red', main = 'Am−Cs−Co spectra', xlab = 'ADC channel [#]', ylab = "Log Counts [#]",cex.main = 1.5)

# rjags model as .bugs file
# peak describe a gaussian
# with an exponential background exp(m(x - mu) + offset)
exp.background.model.4 <- "model{
    #data likelihood
    for (i in 1:length(x)){

         y[i] ~ dpois(points[i])
         
        points[i] <- N/sqrt(2*3.14*sigma*sigma)*exp((-(x[i]-mu)*(x[i]-mu))/(2*sigma*sigma)) + exp(tan(slope.alpha)*(x[i]-mu) + offset)  # line generative model

    }

    #priors for the generative model

    N ~ dunif(40000,70000)            
    mu ~ dunif(4575,4600)        
    sigma ~ dunif(0.1, 6)             
    offset ~ dunif(0,20)        
    slope.alpha ~ dunif(-3.14/2, 0)  


}"

# initialization of the model
peak4.jm <- jags.model(textConnection(exp.background.model.4), 
    data = list(x = peak4$ADC_channel, y = peak4$Counts), n.chain = 3, quiet = TRUE)

# burn-in
update(peak4.jm , 1000, progress.bar="none")

# creates the coda sample
peak4.cs <- coda.samples(model = peak4.jm, variable.names = c("N", "mu", "sigma", "offset", "slope.alpha"), n.iter = 10000, thin = 15, progress.bar="none") 

summary(peak4.cs)

Iterations = 2015:11990
Thinning interval = 15 
Number of chains = 3 
Sample size per chain = 666 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

                  Mean        SD  Naive SE Time-series SE
N            5.572e+04 2.446e+02 5.472e+00      5.473e+00
mu           4.587e+03 1.564e-02 3.498e-04      3.690e-04
offset       2.899e+00 1.793e-02 4.011e-04      4.012e-04
sigma        3.643e+00 1.195e-02 2.673e-04      2.610e-04
slope.alpha -4.865e-03 2.893e-04 6.472e-06      6.618e-06

2. Quantiles for each variable:

                  2.5%        25%        50%        75%      97.5%
N            5.524e+04  5.555e+04  5.571e+04  5.588e+04  5.620e+04
mu           4.587e+03  4.587e+03  4.587e+03  4.587e+03  4.587e+03
offset       2.864e+00  2.887e+00  2.899e+00  2.910e+00  2.934e+00
sigma        3.619e+00  3.636e+00  3.643e+00  3.651e+00  3.667e+00
slope.alpha -5.429e-03 -5.062e-03 -4.866e-03 -4.673e-03 -4.283e-03

# plot of the coda samples
plot(peak4.cs)

# merge the chains on a mcmc objects
peak4_mcmc <- as.mcmc(do.call(rbind, peak4.cs))
#head(peak4_mcmc)

# mean values of the parameters form the mcmc object
peak4_params <- c("N" = mean(peak4_mcmc[,1]),"mu" = mean(peak4_mcmc[,2]), "offset" = mean(peak4_mcmc[,3]), "sigma" = mean(peak4_mcmc[,4]), "alpha" = mean(peak4_mcmc[,5]))
#peak4_params[4]*2.35

# fit of the gaussian peaks + the exponential background
pred_peak4 <-  peak4_params[1] /sqrt(2*3.14*peak4_params[4]*peak4_params[4])*exp((-(peak4$ADC_channel-peak4_params[2])*(peak4$ADC_channel-peak4_params[2]))/(2*peak4_params[4]*peak4_params[4])) +exp(tan(peak4_params[5])*(peak4$ADC_channel - peak4_params[2]) + peak4_params[3]) 

#plot the fit and the measurements
plot(peak4$ADC_channel, peak4$Counts,  col='red', main = 'Am−Cs−Co spectra, Peak 4', xlab = 'ADC channel [#]', ylab = "Counts [#]",cex.main = 1.5)
lines(peak4$ADC_channel, pred_peak4)

COINCIDENCE 60-Co En = 1173.24 keV + 60-Co En = 1332.51 keV

We assume as prior information that the fifth peak is the sum between the two gammas of 60-Co.

This peak will be included in the calibration data set.

# FIFTH PEAK

#range of the fifth peak
peak5 = dataset[!(dataset$ADC_channel<=8300) & (dataset$ADC_channel<= 8700),]
#max4=peak4[which.max(range4$Counts),"ADC_channel"]

# plot the fifth peak
plot(peak5$ADC_channel, peak5$Counts, col='red', main = 'Am−Cs−Co spectra', xlab = 'ADC channel [#]', ylab = "Log Counts [#]",cex.main = 1.5)

# rjags model as .bugs file
# peak describe a gaussian
# with an exponential background exp(m(x - mu) + offset)
exp.background.model.5 <- "model{
    #data likelihood
    for (i in 1:length(x)){

         y[i] ~ dpois(points[i])
         
        points[i] <- N/sqrt(2*3.14*sigma*sigma)*exp((-(x[i]-mu)*(x[i]-mu))/(2*sigma*sigma)) + exp(tan(slope.alpha)*(x[i]-mu) + offset)  # line generative model

    }

    #priors for the generative model

    N ~ dunif(0,10000)            
    mu ~ dunif(8600,8630)        
    sigma ~ dunif(0.1, 10)             
    offset ~ dunif(-5,10)        
    slope.alpha ~ dunif(-3.14/2, 0)  

}"

# initialization of the model
peak5.jm <- jags.model(textConnection(exp.background.model.5), 
    data = list(x = peak5$ADC_channel, y = peak5$Counts), n.chain = 3, quiet = TRUE)
# burn-in
update(peak5.jm , 1000, progress.bar="none")

# coda sample
peak5.cs <- coda.samples(model = peak5.jm, variable.names = c("N", "mu", "sigma", "offset", "slope.alpha"), n.iter = 10000, thin = 15, progress.bar="none") 

summary(peak5.cs)

Iterations = 2015:11990
Thinning interval = 15 
Number of chains = 3 
Sample size per chain = 666 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

                  Mean        SD  Naive SE Time-series SE
N            3.051e+03 5.465e+01 1.223e+00      1.181e+00
mu           8.619e+03 1.250e-01 2.796e-03      2.723e-03
offset      -8.896e-02 8.539e-02 1.910e-03      2.070e-03
sigma        6.676e+00 9.816e-02 2.196e-03      2.121e-03
slope.alpha -1.903e-03 4.352e-04 9.736e-06      1.034e-05

2. Quantiles for each variable:

                  2.5%        25%        50%        75%      97.5%
N            2.942e+03  3.014e+03  3.051e+03  3.088e+03  3.161e+03
mu           8.619e+03  8.619e+03  8.619e+03  8.619e+03  8.619e+03
offset      -2.619e-01 -1.475e-01 -8.707e-02 -2.942e-02  7.308e-02
sigma        6.485e+00  6.609e+00  6.676e+00  6.743e+00  6.861e+00
slope.alpha -2.744e-03 -2.189e-03 -1.898e-03 -1.609e-03 -1.052e-03

# plot of the coda samples
plot(peak5.cs)

# merge the chains in a mcmc object
peak5_mcmc <- as.mcmc(do.call(rbind, peak5.cs))
#head(peak5_mcmc)

# mean values of the parameters
peak5_params <- c("N" = mean(peak5_mcmc[,1]),"mu" = mean(peak5_mcmc[,2]), "offset" = mean(peak5_mcmc[,3]), "sigma" = mean(peak5_mcmc[,4]), "alpha" = mean(peak5_mcmc[,5]))

#  fit of the peak
pred_peak5 <-  peak5_params[1] /sqrt(2*3.14*peak5_params[4]*peak5_params[4])*exp((-(peak5$ADC_channel-peak5_params[2])*(peak5$ADC_channel-peak5_params[2]))/(2*peak5_params[4]*peak5_params[4])) +exp(tan(peak5_params[5])*(peak5$ADC_channel - peak5_params[2]) + peak5_params[3]) 

# plot of data and fit
plot(peak5$ADC_channel, peak5$Counts,  col='red', main = 'Am−Cs−Co spectra, Peak 5', xlab = 'ADC channel [#]', ylab = "Counts [#]",cex.main = 1.5)
lines(peak5$ADC_channel, pred_peak5)

2. Linear Calibration

Assuming a linear response of the detector, as a function of energy, perform a calibration of the detector, associating the centroid of each peak to the nominal value of the detected γ full energy peak

To infer the parameters we use again a JAGS model.

Generative model:

\[\ y = a + b \cdot x\]

• \(a\) = offset

• \(b\) = slope

• \(x\) = centroid the peak

Measurement Model:

\[\text{Norm}(\sigma |E, y ) = \frac{1}{\sqrt{2 \pi}\sigma} \exp \left[ -\frac{(E - y(x))^2}{2\sigma^2}\right] \]

• \(E\): true value of gamma peak, taken from literature

• \(y(x)\): values of the peak from the generative model

• \(sigma\): standard deviation

In this framework we adopt normal prior distributions.

# linear model as .bug file for calibration
# the theoretical value, the fit,  and the measurement will have some discrepancy
# the difference of this two value will follow a normal distribution 
# we also want to infer the value of sigma
linear.model <- "model{
    # Likelihood model for X
    for (i in 1:N){
      E[i] ~ dnorm(f[i], tau) # tau is precision (1 / variance)
          f[i] <- a + b* x[i]
    }
    # Prior model for p
    b ~ dnorm(0, 0.01)
    a ~ dnorm(0, 0.01)
    sigma ~ dunif(0, 100)
    tau <- 1 / (sigma * sigma)

}"

x <- c(2.062e+02 , 2.279e+03, 4.040e+03, 4.587e+03, 8.619e+03) #channel
y <- c(59.54, 661.66, 1173.24, 1332.51, 1173.24 + 1332.51) #true value of energy

# create the jags model
calibration.jags <- jags.model(textConnection(linear.model), 
    data = list(x = x, E = y, N = length(y)), n.chain = 3, quiet =TRUE )

#burn-in
update(calibration.jags , 1000, progress.bar = "none")

# create the coda sample
linear.pred <- coda.samples(model = calibration.jags, variable.names = c("a", "b", "sigma"), n.iter = 10000, thin = 15, progress.bar="none") 

summary(linear.pred)

Iterations = 2015:11990
Thinning interval = 15 
Number of chains = 3 
Sample size per chain = 666 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

         Mean        SD  Naive SE Time-series SE
a     -0.9647 0.9970167 2.231e-02      2.263e-02
b      0.2908 0.0002157 4.826e-06      4.819e-06
sigma  1.0833 0.9256229 2.071e-02      2.598e-02

2. Quantiles for each variable:

         2.5%     25%     50%     75%  97.5%
a     -3.0645 -1.3580 -0.9653 -0.5502 0.9223
b      0.2904  0.2907  0.2908  0.2909 0.2912
sigma  0.3577  0.5720  0.7971  1.2212 3.7649

# plot of the coda samples
plot(linear.pred)

# merge the chain in a mcmc object
calibration_mcmc <- as.mcmc(do.call(rbind, linear.pred))
#head(calibration_mcmc)


# compute the ACF CORRELATIN
calibration.acf <- acf(calibration_mcmc, plot = FALSE) #plot = F will not show the plot of the ACF.

cal_names <- c("b", "a", expression(sigma))

# plot ACF CORRELATION
plot(calibration.acf, type = "s", xlab = "Lag", ylab = NULL,
     ylim = NULL, main = NULL, col = "red", lwd = 2)

#CORRELATION BETWEEN PARAMS
par(mfrow = c(1,3))
for(i in 1:2){
  for(j in (i+1):3){
    points.i <- as.numeric(calibration_mcmc[,i])
points.j <- as.numeric(calibration_mcmc[,j])
plot(points.i, points.j,  col="coral", pch = 1, xlab=cal_names[i], ylab=cal_names[j], cex.lab=1.3, cex.axis=1.3, cex.main=1.3, cex.sub=1.3)
  }
}

# mean values of the calibrations
cal_params <- c(mean(calibration_mcmc[,1]), mean(calibration_mcmc[,2]), mean(calibration_mcmc[,3]))

# grid of channel values
xnew <- 0:10000
n <- length(xnew)

# expectation value of the linear calibration
fit.pred <- cal_params[1] + xnew * cal_params[2]

#extract confidence intervals
confidence_interval_a <- quantile(calibration_mcmc[,"a"], c(0.025, 0.975))
confidence_interval_b <- quantile(calibration_mcmc[,"b"], c(0.025, 0.975))


#plot the expectation value of the linear fit of the calibration
plot(xnew, fit.pred, main = "Calibration of the Ge detector", type = 'l', xlab = "Channel [ch]", ylab = expression("E"* gamma ~ "[keV]"),lwd = 2, col = "dark red")
# measurements
points(x, y, pch = 1, col = "gold", cex = 1, lwd = 2)

# Add legend with formulas
legend_text <- bquote("E"* gamma == .(round(mean(calibration_mcmc[,1]), 2)) + .(round(mean(calibration_mcmc[,2]), 2)) * "* Ch")
sigma_text <- bquote(sigma == .(round(mean(calibration_mcmc[,3]), 2)))
legend("topleft", legend = c(legend_text, sigma_text), bty = "n", cex = 1.2)

# apply the calibration 
dataset$Energy <- dataset$ADC_channel * cal_params[2] + cal_params[1]

#plot of the calibration
plot(dataset$Energy,log(dataset$Counts), type = 's', col='orange', main = 'Calibrated Am−Cs−Co spectra', xlab = expression("E"* gamma ~ "[keV]"), ylab = "Log Counts [#]",cex.main = 1.5)

NA
NA
NA
NA

3. Energy resolution behavior

Using a MCMC method (with either JAGS or stan), study the behaviour of the energy resolution as a function of the photon energy and infer the parameters of eq.1 and 2.

We make an inference using a JAGS model.

Generative model:

\[FWHM = \sqrt{8 ln(2)FwEγ + w_e^2}\]

note that the Fano factor \(F=\dfrac{ \sigma_{sper}^2}{ \sigma_{pois}^2}\). This factor must be introduced since the process giving rise to each individual charge carrier is not independent as the number of ways an atom may be ionized is limited by the discrete electron shells.

Measurement Model:

\[\text{Norm}(\sigma |E, y ) = \frac{1}{\sqrt{2 \pi}\sigma} \exp \left[ -\frac{(E - y(x))^2}{2\sigma^2}\right] \]

• \(E\): true value of gamma peak, taken from literature

• \(y(x)\): values of the peak from the generative model

• \(sigma\): standard deviation

#model as .bug file
#note that w_e in the model should be interpreted as w_e^2 and F as F*w
en.resolution.model <- "model{
    # Likelihood model 
    for (i in 1:N){
      FWHM[i] ~ dnorm(y[i], tau) # tau is precision (1 / variance)
          y[i] <- sqrt(8*log(2)*Fw*E[i] + w_e2) 
    }
    # Prior model for F,w,w_e
    Fw ~ dunif(0, 0.003) 
    #w ~ dunif(0, 0.5)
    w_e2 ~ dunif(0, 10)

    sigma ~ dunif(0, 1)
    tau <- 1 / (sigma * sigma)

}"

# energies of the peak calibrated
x <- c(2.062e+02 , 2.279e+03, 4.040e+03, 4.587e+03)*cal_params[2] + cal_params[1]
#y <- beta_list

# FWHM in channel
y_ch <- c(1.3*2.35,2.379*2.35,3.323*2.35,3.644*2.35) # FWHM best values found

# FWHM in energy
# since it is a range of energies the offset is elided 
y <- y_ch*cal_params[2]


jm <- jags.model(textConnection(en.resolution.model), 
    data = list(E = x, FWHM = y, N = length(y)), n.chain = 3, quiet = TRUE)

#burn-in
update(jm , 1000, progress.bar = "none")

en.resolution <- coda.samples(model = jm, variable.names = c("Fw", "w_e2", "sigma"), n.iter = 10000, thin = 15, progress.bar = "none") 

summary(en.resolution)

Iterations = 2015:11990
Thinning interval = 15 
Number of chains = 3 
Sample size per chain = 666 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

           Mean        SD  Naive SE Time-series SE
Fw    0.0006804 0.0001665 3.725e-06      4.362e-06
sigma 0.2746724 0.1952907 4.369e-03      5.003e-03
w_e2  0.7955765 0.6806435 1.523e-02      1.889e-02

2. Quantiles for each variable:

           2.5%       25%       50%       75%     97.5%
Fw    0.0002572 0.0006169 0.0006976 0.0007597 0.0009701
sigma 0.0752043 0.1301756 0.2077934 0.3561645 0.8081621
w_e2  0.1033196 0.3975490 0.5965618 0.9258164 2.6322145

#pdf("plot_output.pdf")  # Specify the output PDF file name

plot(en.resolution)

#par(mfrow = c(1,3))
#for(i in 1:2){
#  for(j in (i+1):3){
#    points.i <- as.numeric(en.resol_mcmc[,i])
#points.j <- as.numeric(en.resol_mcmc[,j])
#plot(points.i, points.j,  col="coral", pch = 1, xlab=en.resol_names[i], ylab=en.resol_names[j], cex.lab=1.3, cex.axis=1.3, cex.main=1.3, cex.sub=1.3)
#  }
#}

#dev.off()  # Close the PDF device and save the file


en.resol_mcmc <- as.mcmc(do.call(rbind, en.resolution))


#en.resol.acf <- acf(en.resol_mcmc, plot = F) #plot = F will not show the plot of the ACF.

en.resol_names <- c("Fw", expression(sigma), "w_e2")

en.resol_mcmc[,"w_e2"] <- sqrt(en.resol_mcmc[,"w_e2"]) #redefine w_e with its right definition 

#ACF CORRELATION

#plot(en.resol.acf, type = "s", xlab = "Lag", ylab = NULL,
 #    ylim = NULL, main = NULL, col = "red", lwd = 2)

#CORRELATION BETWEEN PARAMS



en.res.params <- c(mean(en.resol_mcmc[,1]), mean(en.resol_mcmc[,2]), mean(en.resol_mcmc[,3]))

#mean(en.resol_mcmc[,1])

E <- seq(from = 0, to = 5000, length.out = 10000)
n <- length(E)

en.resol.fit <- sqrt(8*log(2)*en.res.params[1]*E + en.res.params[3]*en.res.params[3]) 

#plo the data
plot(E, en.resol.fit, main = "Ge detector FWHM vs Energy fit", type = 'l', xlab = expression("E"* gamma ~ "[keV]"), ylab ="FWHM [keV]",lwd = 2, col = "dark red", xlim = c(0,2000))
points(x, y, pch = 1, col = "gold", cex = 2, lwd = 2)

# Add legend with formulas
legend("topleft", legend = c(expression(FWHM == sqrt(8*ln(2)*F*w*E[gamma] + w[e]^2)), bquote(sigma == .(round(mean(en.resol_mcmc[,"sigma"]), 2)))), bty = "n", cex = 1.2)

summary_fit <- data.frame("F" = (mean(en.resol_mcmc[,"Fw"]))/0.003, "w[keV]" = 0.003, "w_e[keV]" = mean(en.resol_mcmc[,"w_e2"]))
names(summary_fit) = c("F", "w[keV]", "w_e[keV]")
summary_fit

NA

Here we have considered \(w\), the electron-hole production energy threshold in germanium, equal to \(3 eV\). To retrieve the Fano factor \(F\) we can just divide \(Fw\) by \(w\). It is difficult to compare rigoursly different Fano factors values with an hypothesis testing since they depend on the Temperature and a lot on the value of \(w\) assigned. However typical values of Fano Factor are smaller than 1.